Ειδικά Μαθήματα Βοτανικής

1 Βασικές προϋποθέσεις

Εννοείται ότι πριν ξεκινήσουμε να κάνουμε το οτιδήποτε, έχουμε δημιουργήσει ένα νέο project στο R-Studio, το όνομα του οποίου - για λόγους τους οποίους θα καταλάβετε αργότερα - ΔΕΝ πρέπει να περιέχει:
1. ελληνικούς χαρακτήρες
2. κενά (αντικαταστήστε τα κενά με _)

Το ίδιο ισχύει ΚΑΙ για το file path του εν λόγω φακέλου (π.χ., όχι ‘E:/Ειδικά Μαθήματα Βοτανικής/SDM,’ αλλά ’E:/Eidika_Mathimata_Botanikis/SDM).

Ένα άλλο κρίσιμο σημείο είναι το εξής: θα χρειαστεί να εγκαταστήσετε τα πακέτα που πρόκειται να χρησιμοποιήσετε σε αυτό το tutorial¹.

2 Εγκατάσταση και φόρτωση βιβλιοθηκών

Όπως κάθε φορά, θα χρειαστεί να εγκαταστήσουμε και να φορτώσουμε τις απαραίτητες βιβλιοθήκες.

## ===========================================================================##
## Install the main packages
## ===========================================================================##
install.packages(c("tidyverse", "rgbif", "raster", "data.table", "rgeos", "rgdal", 
    "gridExtra", "dismo", "biomod2", "sf", "usdm", "speciesgeocodeR", "pacman", "spThin", 
    "CoordinateCleaner", "ggspatial", "scales", "scrubr", "mapr", "rasterVis", "ade4", 
    "viridis", "biogeo"), dependencies = T)
## ===========================================================================##

Ας τις φορτώσουμε:

## ===========================================================================##
## Load them
## ===========================================================================##
pacman::p_load(tidyverse, rgbif, raster, usdm, rgeos, rgdal, usdm, biomod2, gridExtra, 
    countrycode, devtools, sf, spThin, biogeo, speciesgeocodeR, pacman, dismo, CoordinateCleaner, 
    ggspatial, scales, scrubr, mapr, rasterVis, data.table, ade4, viridis)
## ===========================================================================##

3 Διανυσματικά δεδομένα

3.1 Σύνορα χωρών

Πρώτα θα κατεβάσουμε τις γεωγραφικές μεταβλητές για την χώρα που μας ενδιαφέρει (η Ελλάδα στην προκειμένη περίπτωση). Θα πρέπει να εισάγουμε (στην εντολή country στον κώδικα που παρατίθεται πιο κάτω) τον 3-ψήφιο ISO κωδικό για την χώρα μας που τον βρίσκουμε από αυτόν τον ιστότοπο ή με την εντολή ISO_countries <- getData("ISO3") %>% as.data.frame². Εάν θέλουμε να περιορίσουμε την περιοχή μελέτης σε μια μικρότερη περιοχή, τότε θα πρέπει να χρησιμοποιήσουμε ένα άλλο λογισμικό, το Qgis, ώστε να “κόψουμε” την περιοχή που μας ενδιαφέρει. Προς το παρόν, δεν θα χρειαστεί να το κάνουμε αυτό, αλλά θα περιοριστούμε να κατεβάσουμε τα πλήρη δεδομένα για την Ελλάδα: .

## ===========================================================================##
## Download the data for Greece
## ===========================================================================##
## There are three levels available, you are free to explore them
Greece <- getData("GADM", country = "GRC", level = 3)
## ===========================================================================##

Μπορούμε να αναπαραστήσουμε γραφικά τα δεδομένα αυτά.

## ===========================================================================##
## Plot Greece
## ===========================================================================##
plot(Greece)

Είναι αρκετά εύκολο στην προκειμένη περίπτωση να περιορίσουμε την έκταση του πολυγώνου στην Κρήτη, με τις ακόλουθες εντολές:

## ===========================================================================##
## First we need to see what is inside the object we created
## ===========================================================================##
Greece

## class       : SpatialPolygonsDataFrame 
## features    : 326 
## extent      : 19.37236, 29.6457, 34.80069, 41.74801  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 16
## names       : GID_0, NAME_0,   GID_1,                      NAME_1,                                                                                        NL_NAME_1,     GID_2,         NAME_2,                        NL_NAME_2,       GID_3,   NAME_3,          VARNAME_3,                                                     NL_NAME_3, TYPE_3,    ENGTYPE_3, CC_3, ... 
## min values  :   GRC, Greece, GRC.1_1,                      Aegean,                                                                 <U+0386><U+03B8><U+03C9><U+03C2>, GRC.1.1_1,          Athos, <U+0386><U+03B8><U+03C9><U+03C2>, GRC.1.1.1_1,   Abdera, Abelokipi-Menemeni,                              <U+0386><U+03B8><U+03C9><U+03C2>,  Dímos, Municipality,   NA, ... 
## max values  :   GRC, Greece, GRC.8_1, Thessaly and Central Greece, Tessa<U+03BB><U+03AF>a <U+03BA>a<U+03B9> Ste<U+03C1>e<U+03AC> <U+0395><U+03BB><U+03BB><U+03AC>da, GRC.8.2_1, West Macedonia,           Tessa<U+03BB><U+03AF>a, GRC.8.2.9_1, Zografou,           Zografos, Tessa<U+03BB><U+03BF><U+03BD><U+03AF><U+03BA><U+03B7><U+03C2>,  Dímos, Municipality,   NA, ...

names(Greece)

##  [1] "GID_0"     "NAME_0"    "GID_1"     "NAME_1"    "NL_NAME_1" "GID_2"    
##  [7] "NAME_2"    "NL_NAME_2" "GID_3"     "NAME_3"    "VARNAME_3" "NL_NAME_3"
## [13] "TYPE_3"    "ENGTYPE_3" "CC_3"      "HASC_3"

Crete <- subset(Greece, NAME_1 == "Crete")

## ===========================================================================##
## Plot Crete
## ===========================================================================##
plot(Crete)

3.2 GBIF

Στο σημείο αυτό, είμαστε πλέον στην ευχάριστη θέση να μπορούμε να χρησιμοποιήσουμε τα δεδομένα θέσης για τα είδη τα οποία μας ενδιαφέρουν. Υπάρχουν δύο περιπτώσεις, όσον αφορά την απόκτηση τέτοιων δεδομένων:
- να τα έχουμε αποκτήσει πρωτογενώς (δηλαδή να έχουμε εργαστεί στο πεδίο και να έχουμε εντοπίσεις τις θέσεις εμφάνισης του εκάστοτε είδους και των υποπληθυσμών του)
- να τα έχουμε αποκτήσει δευτερογενώς (δηλαδή είτε να έχουμε χρησιμοποιήσει τις διαδικτυακές βάσεις ψηφιοποιημένων - και ενίοτε γεωαναφερμένων - δεδομένων GBIF και BIEN, είτε να έχουμε ψηφιοποιήσει χάρτες εμφάνισης θέσης³)

Το ιδανικό βεβαίως είναι να έχουμε αποκτήσει πρωτογενώς τα στοιχεία αυτά, καθότι στην περίπτωση αυτή, μπορούμε να είμαστε σίγουροι για τις γεωγραφικές συντεταγμένες των υποπληθυσμών του προς μελέτη είδους. Στην πλειονότητα των περιπτώσεων όμως, δεν έχουμε πρόσβαση σε τέτοια δεδομένα, οπότε αναγκαστικά καταφεύγουμε στην δεύτερη λύση και δη στις διαδικτυακές βάσεις δεδομένων. Ευτυχώς για εμάς, υπάρχουν κάποια πακέτα⁴ στην R, όπως το rgbif και το rbien, τα οποία αυτοματοποιούν την διαδικασία μεταφόρτωσης δεδομένων από τις διαδικτυακές αυτές βάσεις⁵.

Εμείς θα εργαστούμε με το πακέτο rgbif προς το παρόν.

Σήμερα – όπως και στην προηγούμενη εργαστηριακή άσκηση – θα ασχοληθούμε ένα εκ των έξι μονοτυπικών γενών της ελληνικής χλωρίδας, την Petromarula pinnata.

Fig 1. The endemic species Petromarula pinnata from Crete.

Πρώτα χρειάζεται να βεβαιωθούμε ότι το προς μελέτη είδος, υπάρχει όντως εντός της διαδικτυακής βάσης την οποία επιθυμούμε να χρησιμοποιήσουμε, προκειμένου να εξάγουμε δεδομένα. Η εντολή name_suggest() ερευνά όλα τα είδη (taxa καλύτερα) τα οποία μοιάζουν με το είδος το οποίο ψάχνουμε και κρατά μόνο το όνομα των taxa το οποίο ταιριάζει ακριβώς με αυτό το οποίο μας ενδιαφέρει⁶.

spp_petromarula <- name_suggest(q = "Petromarula pinnata", rank = "species", limit = 10000)
spp_petromarula

## Records returned [2] 
## No. unique hierarchies [0] 
## Args [q=Petromarula pinnata, limit=10000, rank=species, fields1=key,
##      fields2=canonicalName, fields3=rank] 
## # A tibble: 2 x 3
##        key canonicalName       rank   
##      <int> <chr>               <chr>  
## 1  3167989 Petromarula pinnata SPECIES
## 2 11108006 Petromarula pinnata SPECIES

Στο βήμα αυτό αναγνωρίσαμε τον μοναδικό αριθμό αναφοράς (key) του προς μελέτη είδους και μόνο χρησιμοποιώντας τον αριθμό αυτό μπορούμε να είμαστε βέβαιοι ότι τα δεδομένα θέσης τα οποία θα λάβουμε από την διαδικτυακή αυτή βάση δεδομένων θα αναφέρονται στο είδος το οποίο μας ενδιαφέρει.

Στη συνέχεια, θα δούμε πόσα γεωαναφερμένα δεδομένα θέσης έχουμε στην διάθεση μας.

3.3 Μεταφόρτωση των γεωαναφερμένων δεδομένων θέσης

## ===========================================================================##
## Download data
## ===========================================================================##
sp_occ <- occ_search(taxonKey = spp_petromarula$data$key[1], country = "GR", hasCoordinate = T, 
    limit = 1000)
## ===========================================================================##


## ===========================================================================##
## To prevent any problem with the pathway it is a good practice to remove blank
## spaces from species names
## ===========================================================================##
sp_occ$data$name <- sub(" ", "_", sp_occ$data$name)
## ===========================================================================##


## ===========================================================================##
## Total number of occurrences
## ===========================================================================##
sort(table(sp_occ$data$name), decreasing = T)

## 
## Petromarula_pinnata (L.) A.DC.           Phyteuma_pinnatum L. 
##                            122                              1

Θα πρέπει να τονίσουμε το εξής: Εφ’όσον δεν έχουμε συλλέξει πρωτογενώς τα δεδομένα μας, θα χρειαστεί να ελέγξουμε εάν οι θέσεις εμφάνισης που μεταφορτώσαμε από τις διαδικτυακές βάσεις δεδομένων, συμφωνούν με τα γνωστά όρια εξάπλωσης του είδους που μας ενδιαφέρει. Εάν όντως υπάρχουν κάποια σημεία τα οποία βρίσκονται εκτός της περιοχής εξάπλωσης, θα χρειαστεί να τα αφαιρέσουμε από την ανάλυση.

## ===========================================================================##
## Clean the coordinates
## ===========================================================================##
flags_gbif <- clean_coordinates(x = sp_occ$data, lon = "decimalLongitude", lat = "decimalLatitude", 
    countries = "countryCode", species = "species", tests = c("capitals", "centroids", 
        "equal", "gbif", "zeros", "seas"), seas_ref = Crete)
## ===========================================================================##


## ===========================================================================##
## Keep only those that are OK
## ===========================================================================##
flags_gbif <- flags_gbif %>% mutate(Data = ifelse(.summary == "TRUE", "Clean", "Flagged"), 
    NAME_1 = "Crete", Dataset = "GBIF", longitude = decimalLongitude, latitude = decimalLatitude, 
    scientific_name = "Petromarula pinnata") %>% dplyr::filter(str_detect(Data, "Cle")) %>% 
    dplyr::select(scientific_name, longitude, latitude, Dataset)
## ===========================================================================##


## ===========================================================================##
## Create a duplicate for safekeeping
## ===========================================================================##
flags_gbif_dpl <- flags_gbif
## ===========================================================================##

Ας φτιάξουμε έναν φάκελο ώστε να αποθηκεύσουμε τα δεδομένα μας και ας τα σώσουμε, προκειμένου να μπορέσουμε να τα φορτώσουμε κάποια άλλη στιγμή:

dir.create("./Species", recursive = TRUE, showWarnings = FALSE)
dir.create("./Species/Petromarula pinnata", recursive = TRUE, showWarnings = FALSE)

saveRDS(flags_gbif, "./Species/ Petromarula pinnata.rds")

4 Μεταφόρτωση περιβαλλοντικών μεταβλητών

Τα κλιματικά δεδομένα μας μπορούμε να τα κατεβάσουμε από το WorldClim, ενώ τα υψομετρικά δεδομένα από εδώ και τα δεδομένα σχετικά με την ξηρασία και την εξατμισιοδιαπνοή από εδώ. - Εάν φυσικά δουλεύουμε για μια συγκεκριμένη εργασία (πτυχιακή, μεταπτυχιακή, κτλ) και όχι στα πλαίσια ενός προπτυχιακού μαθήματος επιλογής. Προς το παρόν, θα χρησιμοποιήσουμε κάποιες εντολές για να κατεβάσουμε τα δεδομένα που θέλουμε μέσω της R.

Όπως μπορείτε να διαπιστώσετε, υπάρχουν πάρα πολλά κλιματικά μοντέλα (επί παραδείγματι, CCSM4, HadGEM-2, κτλ) και τέσσερα (4) διαφορετικά κλιματικά σενάρια που έχουν να κάνουν με τη συγκέντρωση των αερίων του θερμοκηπίου στην ατμόσφαιρα εντός του 21^ου αιώνα { ενδεικτική βιβλιογραφία ή άλλη ενδεικτική βιβλιογραφία ή ακόμα μια ενδεικτική βιβλιογραφία ή τελευταία ενδεικτική βιβλιογραφία}⁷.

Συνεπώς, δεν μπορούμε να χρησιμοποιήσουμε μόνο ένα κλιματικό μοντέλο, πόσω μάλλον ένα κλιματικό σενάριο. Στο σημείο αυτό, είναι εύλογο να αρχίσετε να πανικοβάλεστε⁸. Ευτυχώς όμως, δεν χρειάζεται και ούτε είναι επιστημονικά σωστό, να χρησιμοποιήσουμε όλα τα κλιματικά μοντέλα. Αναλόγως της περιοχής στην οποία εργαζόμαστε, κάποια κλιματικά μοντέλα αναπαριστούν καλύτερα τις περιβαλλοντικές συνθήκες σε σχέση με κάποια άλλα (McSweeney et al., 2015). Καθώς θα διαπιστώσετε (όταν με το καλό διαβάσετε την προαναφερθείσα εργασία), ακόμα και έτσι, είναι πολλά τα κλιματικά μοντέλα για την περιοχή μας. Όμως, μόνο 3 από αυτά έχουν δεδομένα και για τα 4 κλιματικά σενάρια (BCC, HadGEM, CCSM4). Εμείς χάριν συντομίας, θα εργαστούμε όμως μόνο με ένα (1) κλιματικό μοντέλο και ενα (1) κλιματικό σενάριο.

## ============================================================================##
## First, let's download current climate data for the whole planet
## ============================================================================##
current_climate <- getData("worldclim", var = "bio", res = 2.5)
## ============================================================================##


## ============================================================================##
## Then, let's download future climate data
## ============================================================================##
future_climate <- getData("CMIP5", var = "bio", res = 2.5, rcp = 85, model = "CC", 
    year = 70)

## ============================================================================##

Ας αλλάξουμε τα ονόματα των περιβαλλοντικών μεταβλητών.

names(current_climate) <- c("bio_1", "bio_2", "bio_3", "bio_4", "bio_5", "bio_6", 
    "bio_7", "bio_8", "bio_9", "bio_10", "bio_11", "bio_12", "bio_13", "bio_14", 
    "bio_15", "bio_16", "bio_17", "bio_18", "bio_19")

names(future_climate) <- c("bio_1", "bio_2", "bio_3", "bio_4", "bio_5", "bio_6", 
    "bio_7", "bio_8", "bio_9", "bio_10", "bio_11", "bio_12", "bio_13", "bio_14", 
    "bio_15", "bio_16", "bio_17", "bio_18", "bio_19")

Καθώς τα κλιματικά μας δεδομένα όπως διαπιστώσατε έχουν εξαιρετικά μεγάλο εύρος τιμών, θα χρειαστεί να τα μετατρέψουμε σε ^oC όσον αφορά την θερμοκρασία (bio_1). Εδώ μπορείτε να διαπιστώσετε το γιατί⁹. Αυτό γίνεται με την ακόλουθη εντολή:

gain(current_climate$bio_1) = 0.1

gain(future_climate$bio_1) = 0.1

Στη συνέχεια, θα κατεβάσουμε τα υψομετρικά δεδομένα για την χώρα που μας ενδιαφέρει:

altitude <- getData("alt", country = "GRC", mask = TRUE)

4.1 Κοπή των περιβαλλοντικών δεδομένων στην έκταση της Ελλάδας

Στη συνέχεια, θα ‘κόψουμε’ τις περιβαλλοντικές μας μεταβλητές στο σχήμα της Ελλάδας και ακολούθως θα δημιουργήσουμε ένα νέο αρχείο που θα περιέχει τις μεταβλητές αυτές και την περιοχή μελέτης μας:

## ============================================================================##
## First for Greece
## ============================================================================##
current_climate_crop <- crop(current_climate, Greece, snap = "in") %>% mask(., Greece)

altitude <- crop(altitude, Greece, snap = "in") %>% mask(., Greece) %>% 
## We are forcing the extent of altitude to be exactly the same as that of the
## climate data
resample(., current_climate_crop, "bilinear")

future_climate <- crop(future_climate, Greece, snap = "in") %>% mask(., Greece)

## ============================================================================##


## ============================================================================##
## Merge (stack) climate and altitude
## ============================================================================##
current_climate_crop <- stack(current_climate_crop, altitude)
future_climate_crop <- stack(future_climate, altitude)
## ============================================================================##


## ============================================================================##
## Then for Crete
## ============================================================================##
current_climate_Crete <- crop(current_climate_crop, Crete, snap = "in") %>% mask(., 
    Crete)

future_climate_Crete <- crop(future_climate_crop, Crete, snap = "in") %>% mask(., 
    Crete)

## ============================================================================##

4.2 Γραφική απεικόνιση

Μπορούμε να αναπαραστήσουμε γραφικά τα αποτελέσματα ως εξής:

gplot(current_climate_Crete[[1]]) + geom_raster(aes(fill = value)) + scale_fill_gradientn(colours = c("brown", 
    "red", "yellow", "darkgreen", "green")) + coord_equal() + xlab("Longitude") + 
    ylab("Latitude") + ggtitle("Temperature")

histogram(current_climate_Crete[[1]])

levelplot(current_climate_Crete[[1]])

Κάντε το ίδιο για την Αττική, τον Άθω και την Θεσσαλία.

4.3 Βασικά στατιστικά στοιχεία

Στη συνέχεια, θα υπολογίσουμε κάποια βασικά στατιστικά στοιχεία για κάθε αρχείο το οποίο δημιουργήσαμε:

cellStats(current_climate_crop, "quantile")

##      bio_1 bio_2 bio_3   bio_4 bio_5 bio_6 bio_7 bio_8 bio_9 bio_10 bio_11
## 0%     2.4    55    27 4730.00   180   -95   182   -37    94    106    -53
## 25%   12.2    92    33 5988.75   280   -11   262    60   205    207     37
## 50%   14.5   101    34 6516.00   298    12   288    81   229    230     60
##      bio_12 bio_13 bio_14 bio_15 bio_16 bio_17 bio_18 bio_19 GRC_msk_alt
## 0%      385  52.00      0     20  138.0      2      2    110   -2.484375
## 25%     542  74.00      8     32  206.0     34     39    184  111.520996
## 50%     653 103.00     15     50  269.5     65     69    250  289.577637
##  [ reached getOption("max.print") -- omitted 2 rows ]

cellStats(current_climate_crop, "mean")

##       bio_1       bio_2       bio_3       bio_4       bio_5       bio_6 
##    14.03543    97.27204    34.09507  6455.04935   295.72570    15.26490 
##       bio_7       bio_8       bio_9      bio_10      bio_11      bio_12 
##   280.46080    80.48347   222.14066   224.08342    60.01553   695.86060 
##      bio_13      bio_14      bio_15      bio_16      bio_17      bio_18 
##   110.14041    15.94489    51.10471   296.70353    61.12325    66.01879 
##      bio_19 GRC_msk_alt 
##   269.81864   430.21763

cellStats(current_climate_crop, "sd")

##       bio_1       bio_2       bio_3       bio_4       bio_5       bio_6 
##    2.854903   14.953890    2.226295  671.676911   25.702217   37.634237 
##       bio_7       bio_8       bio_9      bio_10      bio_11      bio_12 
##   35.202266   31.311853   27.093454   26.693768   33.891841  189.062736 
##      bio_13      bio_14      bio_15      bio_16      bio_17      bio_18 
##   41.171401    9.894798   20.722858  109.157622   32.687959   35.129005 
##      bio_19 GRC_msk_alt 
##  102.446221  401.097651

Κάντε το ίδιο για την Αττική, τον Άθω και την Θεσσαλία.

Μπορούμε να δημιουργήσουμε υποσύνολα των περιοχών που μας ενδιαφέρουν, ώστε να ανταποκρίνονται σε συγκεκριμένες κλιματικές (ή άλλες συνθήκες):

subset_GR <- current_climate_crop[[1]] < 15 & current_climate_crop[[1]] > 10.1

plot(subset_GR)

Δοκιμάστε το και για άλλες περιοχές και μεταβλητές.

4.4 Σύγκριση

Τώρα, θα πρέπει να ελέγξουμε εάν τα υψομετρικά δεδομένα έχουν την ίδια ανάλυση (resolution), μέγεθος (nrow, ncol) και έκταση (bbox), με τα περιβαλλοντικά μας δεδομένα:

altitude

## class      : RasterLayer 
## dimensions : 165, 246, 40590  (nrow, ncol, ncell)
## resolution : 0.04166667, 0.04166667  (x, y)
## extent     : 19.375, 29.625, 34.83333, 41.70833  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : GRC_msk_alt 
## values     : -2.484375, 2346.14  (min, max)

current_climate_crop

## class      : RasterStack 
## dimensions : 165, 246, 40590, 20  (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667  (x, y)
## extent     : 19.375, 29.625, 34.83333, 41.70833  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      :       bio_1,       bio_2,       bio_3,       bio_4,       bio_5,       bio_6,       bio_7,       bio_8,       bio_9,      bio_10,      bio_11,      bio_12,      bio_13,      bio_14,      bio_15, ... 
## min values :    2.400000,   55.000000,   27.000000, 4730.000000,  180.000000,  -95.000000,  182.000000,  -37.000000,   94.000000,  106.000000,  -53.000000,  385.000000,   52.000000,    0.000000,   20.000000, ... 
## max values :       19.30,      123.00,       41.00,     7606.00,      348.00,      104.00,      341.00,      155.00,      270.00,      270.00,      133.00,     1336.00,      226.00,       49.00,      102.00, ...

4.5 Resolution change

Στην περίπτωση κατά την οποία τα δύο σετ δεδομένων έχουν διαφορετική ανάλυση, έκταση και μέγεθος, θα χρειαστεί να αλλάξουμε την ανάλυση των υψομετρικών δεδομένων¹⁰.

agg_altitude <- aggregate(altitude, 
                 5, ## Why choose 5 and not any other number?
                 fun = mean)

4.6 Pixel size

Μετά, θα αλλάξουμε το μέγεθος των pixel των υψομετρικών δεδομένων για να συμφωνεί με αυτό των περιβαλλοντικών δεδομένων:

altitude_rsmp <- resample(agg_altitude, current_climate_crop, method = "bilinear")

4.7 Final stack

Τέλος, θα τα ενώσουμε σε ένα αρχείο, τόσο για το παρόν, όσο και για το μέλλον

## ============================================================================##
## Merge (stack) climate and altitude
## ============================================================================##
current_climate_crop <- stack(current_climate_crop, altitude_rsmp)
future_climate_crop <- stack(future_climate_crop, altitude_rsmp)
## ============================================================================##

4.8 Αποθήκευση σε συμπιεσμένη μορφή

Καλό είναι να αποθηκεύσετε τα αρχεία αυτά υπό την μορφή Rds και να τα αφαιρέσετε από την μνήμη της R11. Αυτό το κάνουμε, προκειμένου να μπορούμε εύκολα να το φορτώνουμε στην R και να μην υπάρχουν τα δεδομένα αυτά συνεχώς στην μνήμη της session που έχουμε ανοικτή, διότι αυξάνουν αρκετά τον χρόνο απόκρισης της (καθυστερεί πολύ η R). Μπορείτε ανά πάσα στιγμή να τα φορτώσετε στην R, με την ακόλουθη εντολή:

## Creating a folder to save our data
dir.create("./Data/", recursive = TRUE, showWarnings = FALSE)

## Save our data in rds format
saveRDS(current_climate_crop, "./Data/Environmental variables and Elevation.rds")
saveRDS(future_climate_crop, "./Data/Environmental variables and Elevation 2070.rds")

## Load our data from rds format
var_GR <- readRDS("./Data/Environmental variables and Elevation.rds")
var_GR_70 <- readRDS("./Data/Environmental variables and Elevation 2070.rds")

4.9 Έλεγχος συγγραμικότητας και Variance Inflation Factor (VIF)

Όπως θα γνωρίζετε, ΔΕΝ μπορούμε να χρησιμοποιήσουμε σε μια ανάλυση δύο μεταβλητές, οι οποίες εμφανίζουν συντελεστή συσχέτισης \(> 0.7\)¹¹, καθότι στην προκειμένη περίπτωση οι δύο αυτές μεταβλητές είναι συγγραμικές και ως εκ τούτου εάν τις συμπεριλάβουμε και τις δύο στην ανάλυση μας, τότε δεν θα γνωρίζουμε σε ποιά εκ των δύο θα οφείλεται το αποτέλεσμα της ανάλυσης μας.

Προκειμένου να ξεπεράσουμε αυτόν τον σκόπελο, θα τρέξουμε την ακόλουθη εντολή¹²:
.

uncorrelated_vars <- vifcor(current_climate_Crete %>% as.data.frame(), th = 0.7)
uncorrelated_vars <- uncorrelated_vars@results$Variables[uncorrelated_vars@results$VIF <= 
    5]

current_climate_Crete <- stack(subset(current_climate_Crete, uncorrelated_vars))

future_climate_Crete <- stack(subset(future_climate_Crete, uncorrelated_vars))

#Δεδομένα θέσης και περιβαλλοντικές μεταβλητές

Πλέον μπορούμε να συνδέσουμε τις γεωαναφερμένες θέσεις εμφάνισης του είδους που μας ενδιαφέρει με δεδομένα ψηφιδοπλέγματος. Στο σημείο αυτό, μπορούμε να εξάγουμε τα περιβαλλοντικά δεδομένα για κάθε μια εκ των θέσεων εμφάνισης του είδους που μας ενδιαφέρει με μια μόνο εντολή, αναδεικνύοντας κατ’αυτόν τον τρόπο το πλεονέκτημα του να έχουμε ενώσει προηγουμένως τα raster αρχεία μας σε ένα RasterStack.

flags_gbif <- flags_gbif %>% as.data.frame()

coordinates(flags_gbif) <- c("longitude", "latitude")

pts.clim <- raster::extract(current_climate_Crete, flags_gbif, method = "bilinear")

sp_occ_clim <- data.frame(cbind(coordinates(flags_gbif), pts.clim, flags_gbif@data))

Και τώρα, ήρθε η στιγμή να οπτικοποιήσουμε τα αποτελέσματα μας και δη, την παρουσία του είδους μας στην περιοχή μελέτης:

map.ext <- extent(Crete)

plot(current_climate_crop[[20]], col = terrain.colors(100), ext = map.ext)

plot(flags_gbif, pch = 16, cex = 0.5, add = T)

plot(Crete, add = T)

Στο σημείο αυτό, θα ήθελα να σας τονίσω ότι ορισμένες φορές, πέρα από τους ελέγχους που διενεργήσαμε έως τώρα, ίσως χρειαστεί να προχωρήσουμε σε μια αραίωση (thinning) των θέσεων εμφάνισης. Αυτό μπορεί να είναι απαραίτητο όταν έχουμε πάρα πολλά σημεία τα οποία είναι αρκετά κοντά το ένα στο άλλο¹³. Σε κάθε περίπτωση, το επιστημονικά ορθό είναι να υπάρχει μια (1) θέση εμφάνισης ανά km². Ένα πακέτο που μπορεί να μας βοηθήσει είναι το spThin¹⁴.

## Min. distance: 1000 m Number of reps: 1000
spoccs <- flags_gbif_dpl %>% as.data.frame() %>% dplyr::rename(x = longitude, y = latitude, 
    Species = scientific_name) %>% addmainfields(species = "Species")

thinned <- thin(loc.data = spoccs, lat.col = "y", long.col = "x", spec.col = "Species", 
    thin.par = 5, reps = 1000, verbose = T, out.dir = "./Thinned", locs.thinned.list.return = TRUE, 
    write.files = TRUE, out.base = "Petromarula pinnata thinned")

## ********************************************** 
##  Beginning Spatial Thinning.
##  Script Started at: Sun Jun 06 18:46:29 2021
## lat.long.thin.count
##   52 
## 1000 
## [1] "Maximum number of records after thinning: 52"
## [1] "Number of data.frames with max records: 1000"
## [1] "Writing new *.csv files"
## [1] "Writing file: ./Thinned/Petromarula pinnata thinned_thin1.csv"
## [1] "Writing file: ./Thinned/Petromarula pinnata thinned_thin2.csv"
## [1] "Writing file: ./Thinned/Petromarula pinnata thinned_thin3.csv"
## [1] "Writing file: ./Thinned/Petromarula pinnata thinned_thin4.csv"
## [1] "Writing file: ./Thinned/Petromarula pinnata thinned_thin5.csv"

5 Μοντελοποίηση

Όλα τα απαιτούμενα δεδομένα και αρχεία θα τα βρείτε στον ιστότοπο του μαθήματος στο εδώ

5.1 Εισαγωγή

Είμαστε πλεόν στην ευχάριστη θέση να έχουμε στην διάθεση μας όλα τα δεδομένα τα οποία είναι απαραίτητα για να προχωρήσουμε στο κυρίως σκέλος της ανάλυσης μας. Κανονικά θα έπρεπε να χρησιμοποιήσουμε πολλά μοντέλα, τα οποία ανήκουν στην ίδια ομάδα (Παρουσίας/Απουσίας ή Ψεύδο-απουσίας - Presence/Absence: P/A), αλλά και σε άλλη ομάδα (Μόνο Παρουσίας - Presence Only: P/O). Κάθε ένα εκ των μοντέλων αυτών έχει τα πλεονεκτήματα και τα μειονεκτήματα του¹⁵. Στο σημείο αυτό, πρέπει να τονίσουμε ότι είναι εντελώς λάθος να χρησιμοποιούμε δύο μοντέλα¹⁶ τα οποία ανήκουν σε διαφορετικές ομάδες (εν αντιθέσει με ότι προτείνεται στο παράδειγμα του προαναφερθέντος βιβλίου), αλλά να χρησιμοποιεί πολλά μοντέλα τα οποία ανήκουν στην ίδια ομάδα (π.χ., P/A), προκειμένου να μπορεί να βγάλει ασφαλή συμπεράσματα σχετικά με το ερώτημα το οποίο καλείται να απάντησει και να είναι σίγουρος ότι δεν πρόκειται περί στατιστικών κατασκευασμάτων (Araújo & New, 2007).

5.2 Φόρτωση των δεδομένων θέσης

Τώρα θα φορτώσουμε ένα από τα set δεδομένων¹⁷ τα οποία έχουν προέλθει από το thinning των θέσεων εμφάνισης¹⁸.

thinned_1 <- readxl::read_excel("G:/Academic/Lessons/EMB/Labs/Labs/Fourth_Lab/Fourth_Lab/Thinned/Petromarula pinnata thinned_thin1.xlsx")

5.3 Μορφοποιήση των δεδομένων

Το πρώτο πράγμα το οποίο θα χρειαστεί να κάνουμε, είναι να μορφοποιήσουμε με τον κατάλληλο τρόπο τα δεδομένα μας, ώστε να μπορέσουμε να τρέξουμε τις εντολές από τα πακέτα biomod και ecospat. Το δεύτερο πακέτο θα το χρησιμοποιήσουμε στην περίπτωση όπου οι θέσεις εμφάνισης του είδους που μας ενδιαφέρει είναι \(<50\)¹⁹.

Ένα άλλο σημείο το οποιο πρέπει να προσέξουμε, έχει να κάνει με τον αριθμό των απουσιών (absences - στην πλειονότητα των περιπτώσεων βέβαια πρόκειται για ψευδοαπουσίες -pseudoabsences, καθώς στην πραγματικότητα δεν έχουμε δεδομένα απουσίας του εκάστοτε είδους). Ο αριθμός αυτός πρέπει να είναι ίσος με τον αριθμό των παρουσιών²⁰.

Τέλος, ένα τρίτο, κρίσιμο σημείο είναι ο αριθμός των set των pseudoabsences που θα περιλαμβάνει το μοντέλο μας, προκειμένου να αποφύγουμε τον σκόπελο της μεροληπτικής δειγματοληψίας (sampling bias)²¹.

sp_formated_data <- BIOMOD_FormatingData(resp.var = rep(1, nrow( thinned_1 ) ),
                                         
                                    ## Our rasterStack with the uncorrelated
                                    ## variables
                                    expl.var = current_climate_Crete,
                                    
                                    ## The names of our coords
                                    resp.xy = thinned_1[,c('x', 'y')],
                                    
                                    ## Species' name
                                    resp.name = "Petromarula pinnata",
                                    
                                    ## Number of pseudoabsences' sets
                                    PA.nb.rep = 2,
                                    
                                    ## Number equaling that of the presences
                                    PA.nb.absences = nrow(thinned_1), 
                                    
                                    ## You can change this to 'disk' if you'd
                                    ## like
                                    PA.strategy = 'random',
                                    
                                    ## Min distance from occurrences
                                    PA.dist.min  =  NULL,
                                    
                                    ## Max distance from occurrences
                                    PA.dist.max  =  NULL)

## 
## -=-=-=-=-=-=-=-=-=-=-= Petromarula pinnata Data Formating -=-=-=-=-=-=-=-=-=-=-=
## 
##  Response variable name was converted into Petromarula.pinnata
##       ! No data has been set aside for modeling evaluation
##    > Pseudo Absences Selection checkings...
##    > random pseudo absences selection
##    > Pseudo absences are selected in explanatory variables
##          ! Some NAs have been automatically removed from your data
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Done -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

sp_formated_data

## 
## -=-=-=-=-=-=-=-=-=-=-=-=-= 'BIOMOD.formated.data.PA' -=-=-=-=-=-=-=-=-=-=-=-=-=
## 
## sp.name =  Petromarula.pinnata
## 
##   31 presences,  0 true absences and  75 undifined points in dataset
## 
## 
##   4 explanatory variables
## 
##      bio_3           bio_5           bio_13          bio_15     
##  Min.   :29.00   Min.   :222.0   Min.   :111.0   Min.   :80.00  
##  1st Qu.:30.00   1st Qu.:269.0   1st Qu.:138.2   1st Qu.:83.00  
##  Median :31.00   Median :277.0   Median :153.0   Median :84.00  
##  Mean   :31.24   Mean   :275.0   Mean   :160.3   Mean   :84.02  
##  3rd Qu.:32.00   3rd Qu.:286.8   3rd Qu.:184.8   3rd Qu.:85.00  
##  Max.   :33.00   Max.   :301.0   Max.   :213.0   Max.   :89.00  
## 
## 
##  2 Pseudo Absences dataset available ( PA1 PA2 ) with  39 
## absences in each (true abs + pseudo abs)
## 
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

Ας αναπαραστήσουμε γραφικά ότι φτιάξαμε στο προηγούμενο βήμα:

plot(sp_formated_data)

Η εικόνα αυτή μας δείχνει που βρίσκονται οι ψευδοαπουσίες (pseudoabsences) στο set των ψευδοαπουσιών σε σχέση με τις θέσεις εμφάνισης. Όπως αναμενόταν, οι ψευδοαπουσίες είναι τυχαία κατανεμημένες.

5.4 Παραμετροποιήση των μοντέλων

Πλέον μπορούμε να παραμετροποιήσουμε τα μοντέλα που θέλουμε να τρέξουμε:

modelling_options <- BIOMOD_ModelingOptions()

Τώρα, προσδιορίζουμε ποιά μοντέλα θα τρέξουμε, καθώς και πολλές άλλες λεπτομέρειες:

## Always change the name of this object, so as to match the name of your species
Petromarula_pinnata_models <- BIOMOD_Modeling( data = sp_formated_data,
                                       
                                       
                                       ## Here, we indicate which models we are
                                       ## going to run (remember that they need
                                       ## to be from the same group P/A ή P/O).
                                       ## Here, we are going to use some models
                                       ## belonging to the P/A group
                                       ## ('GLM','GBM','GAM','CTA','ANN','SRE',
                                       ## 'FDA','MARS','RF')
                                       models  =  c('RF', 'CTA'),
                                       
                                       models.options = modelling_options,
                                       
                                       ## How many times we will run the
                                       ## analysis in order to evaluate its'
                                       ## results [min: 10 - usually: 100]
                                       NbRunEval = 1,
                                       
                                       ## Training-Testing data split
                                       DataSplit = 80, 
                                       
                                       ## How many times we will run the
                                       ## analysis in order to evaluate the
                                       ## importance of each IV [min: 10]
                                       VarImport = 1,
                                       
                                       ## Evaluation criteria
                                       models.eval.meth  =  c('TSS','ROC','KAPPA'), 
                                       
                                       ## Whether or not we will save the results
                                       SaveObj  =  TRUE,
                                       
                                       ## Should the models be re-scaled in 0-1000 range?
                                       rescal.all.models  =  TRUE,
                                       
                                       ## Should our models be weighted based on all the data?
                                       do.full.models  =  TRUE, 
                                       
                                       ## The name of our analysis
                                       modeling.id  =  "allmodels")

Ας φτιάξουμε τους φακέλους όπου θα αποθηκεύσουμε τα μοντέλα μας:

## Create some folders to store our results
dir.create("./BIOMOD2 Modelling Outpout/", recursive = TRUE, showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/", recursive = TRUE, showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata", recursive = TRUE, 
    showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/ Petromarula_pinnata/Current", recursive = TRUE, 
    showWarnings = FALSE)


## Save the results in rds format
saveRDS(Petromarula_pinnata_models, "./BIOMOD2 Modelling Outpout/Species/ Petromarula_pinnata/Current/ Petromarula_pinnata models.rds")

5.5 Αξιολόγηση των αρχικών μοντέλων

Ας δημιουργήσουμε ένα αντικείμενο το οποίο θα περιέχει την αξιολόγηση των μοντέλων:

Petromarula_pinnata_scores <- get_evaluations(Petromarula_pinnata_models)

5.6 Σημαντικότητα μεταβλητών

Ας δούμε τη σημαντικότητα των μεταβλητών για κάθε μοντέλο (RF, CTA), κάθε σετ δεδομένων (PA1-2) και κάθε επανάληψη (Run, Full):

(Petromarula_pinnata_var_import <- get_variables_importance(Petromarula_pinnata_models))

Ας δούμε τη σημαντικότητα των μεταβλητών για κάθε μοντέλο (RF, CTA) συνολικά (δηλαδή τον μέσο όρο):

options(digits = 2)
apply(Petromarula_pinnata_var_import, c(1, 2), mean)

##          RF   CTA
## bio_3  0.17 0.355
## bio_5  0.33 0.431
## bio_13 0.37 0.632
## bio_15 0.17 0.086

5.7 Ensemble modelling

Τώρα θα τρέξουμε την ανάλυση που θα ενώνει τα αποτελέσματα των δύο μοντέλων (ονομάζεται ensemble modelling αυτή η τεχνική ανάλυσης):

Petromarula_pinnata_ensemble_models <- BIOMOD_EnsembleModeling( modeling.output = Petromarula_pinnata_models,
                                                   em.by = 'all',
                                                   chosen.models  =  'all',
                                                   
                                                   # You can change that
                                                   eval.metric = 'TSS', 
                                                   
                                                   # Keep only the models with
                                                   # TSS >= 0.8
                                                   eval.metric.quality.threshold = 0.8, 
                                                   models.eval.meth = c('KAPPA','TSS','ROC'),
                                                   prob.mean  =  T, 
                                                   committee.averaging = TRUE,
                                                   prob.mean.weight  =  T, 
                                                   prob.mean.weight.decay  =  'proportional',
                                                   VarImport = 0 )

dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Ensemble", 
    recursive = TRUE, showWarnings = FALSE)

saveRDS(Petromarula_pinnata_ensemble_models, "./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Ensemble/Petromarula_pinnata ensemble.rds")

5.8 Αξιολόγηση των συνολικών μοντέλων (ensemble modelling)

Ας δημιουργήσουμε ένα αντικείμενο το οποίο θα περιέχει την αξιολόγηση των μοντέλων από το ensemble modelling. Όσο πιο κοντά στη μονάδα είναι οι τιμές για τους δείκτες TSS, ROC και KAPPA, τόσο καλύτερα²²:

options(digits = 4)
(Petromarula_pinnata_ensemble_models_scores <- get_evaluations(Petromarula_pinnata_ensemble_models))

## $Petromarula.pinnata_EMmeanByTSS_mergedAlgo_mergedRun_mergedData
##       Testing.data Cutoff Sensitivity Specificity
## KAPPA        0.911  616.0       96.77          96
## TSS          0.928  616.0       96.77          96
## ROC          0.994  618.5       96.77          96
## 
## $Petromarula.pinnata_EMcaByTSS_mergedAlgo_mergedRun_mergedData
##       Testing.data Cutoff Sensitivity Specificity
## KAPPA        0.909    747       93.55       97.33
## TSS          0.909    747       93.55       97.33
## ROC          0.969    750       93.55       97.33
## 
## $Petromarula.pinnata_EMwmeanByTSS_mergedAlgo_mergedRun_mergedData
##       Testing.data Cutoff Sensitivity Specificity
## KAPPA        0.911  611.0       96.77          96
## TSS          0.928  611.0       96.77          96
## ROC          0.994  611.5       96.77          96

## ===========================================================## Vizualisation
## ===========================================================##

g2 <- models_scores_graph(Petromarula_pinnata_ensemble_models, by = "models", metrics = c("TSS", 
    "KAPPA"))

5.9 Προβολή στο μέλλον

Ας τρέξουμε την ανάλυση για το πώς θα αλλάξει η κατάσταση του είδους που μας ενδιαφέρει στο μέλλον. Αρχικά θα τρέξουμε την ανάλυση για κάθε μοντέλο ξεχωριστά (BIOMOD_Projection) και συνολικά (BIOMOD_EnsembleForecasting) για τις παρούσες συνθήκες και στη συνέχεια για τις μελλοντικές συνθήκες (με τις ίδιες εντολές, αλλά θα αλλάξουμε το raster αρχείο στην εντολή new.env):

Petromarula_pinnata_models_proj_current <- BIOMOD_Projection(modeling.output = Petromarula_pinnata_models, 
    new.env = current_climate_Crete, proj.name = "current", binary.meth = "TSS", 
    selected.models = "all", build.clamping.mask = TRUE, output.format = ".img", 
    do.stack = T)

Petromarula_pinnata_ensemble_models_proj_current <- BIOMOD_EnsembleForecasting(EM.output = Petromarula_pinnata_ensemble_models, 
    projection.output = Petromarula_pinnata_models_proj_current, binary.meth = "TSS", 
    output.format = ".img", do.stack = T)

dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/", 
    recursive = TRUE, showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Individual", 
    recursive = TRUE, showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Ensemble", 
    recursive = TRUE, showWarnings = FALSE)

saveRDS(Petromarula_pinnata_models_proj_current, "./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Individual/Petromarula_pinnata Current Individual Projections.rds")

saveRDS(Petromarula_pinnata_ensemble_models_proj_current, "./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Ensemble/Petromarula_pinnata Current Ensemble Projections.rds")

Petromarula_pinnata_proj_2070_cc85 <- BIOMOD_Projection(modeling.output = Petromarula_pinnata_models, 
    new.env = future_climate_Crete, proj.name = "2070_CC85", binary.meth = "TSS", 
    output.format = ".img", do.stack = T)

dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070", recursive = TRUE, 
    showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/", 
    recursive = TRUE, showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Individual", 
    recursive = TRUE, showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Ensemble", 
    recursive = TRUE, showWarnings = FALSE)

saveRDS(Petromarula_pinnata_proj_2070_cc85, "./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Individual/Petromarula_pinnata 2070 Individual Projections.rds")

Petromarula_pinnata_ensemble_models_proj_2070_cc85 <- BIOMOD_EnsembleForecasting(EM.output = Petromarula_pinnata_ensemble_models, 
    projection.output = Petromarula_pinnata_proj_2070_cc85, binary.meth = "TSS", 
    output.format = ".img", do.stack = T)

saveRDS(Petromarula_pinnata_ensemble_models_proj_2070_cc85, "./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Ensemble/Petromarula_pinnata 2070 Ensemble Projections.rds")

5.10 Οπτικοποιήση των συνολικών μοντέλων

Τέλος, ας οπτικοποιήσουμε το αποτέλεσμα της ανάλυσης μας. Πρώτα θα θα δημιουργήσουμε ένα αρχείο το οποίο θα περιέχει την αξιολόγηση της κατάστασης της Petromarula pinnata μέσω της τεχνικής ensemble modelling. Στη συνέχεια, θα κρατήσουμε μόνο τα raster layers για το mean και το weighted mean. Τέλος, θα δημιουργήσουμε ένα αρχείο το οποίο θα περιέχει μόνο τις συντεταγμένες των θέσεων εμφάνισης του είδους που μας ενδιαφέρει:

stk_Petromarula_pinnata_ef_2070_cc85 <- get_predictions(Petromarula_pinnata_ensemble_models_proj_2070_cc85)

stk_Petromarula_pinnata_ef_2070_cc85 <- subset(stk_Petromarula_pinnata_ef_2070_cc85, 
    grep("EMmean|EMwmean", names(stk_Petromarula_pinnata_ef_2070_cc85)))

names(stk_Petromarula_pinnata_ef_2070_cc85) <- sapply(strsplit(names(stk_Petromarula_pinnata_ef_2070_cc85), 
    "_"), getElement, 2)

Petromarula_pinnata_coords <- thinned_1 %>% dplyr::select(x, y)

Στη συνέχεια θα οπτικοποιήσουμε το αποτέλεσμα της ανάλυσης μας έως αυτό το σημείο:

library(viridis)
levelplot(stk_Petromarula_pinnata_ef_2070_cc85, main = "Petromarula pinnata ensemble projections\nin 2070 with cc85", 
    col.regions = viridis(256)) + layer(sp.points(SpatialPoints(Petromarula_pinnata_coords), 
    pch = 20, col = 1))

5.11 Μεταβολή περιοχής εξάπλωσης

Έχουμε φτάσει πλέον στο τελικό σημείο της ανάλυσης μας, όπου θα φορτώσουμε τα αποτελέσματα της ανάλυσης για το mean και το weighted mean της Petromarula pinnata και στη συνέχεια θα υπολογίσουμε και θα οπτικοποιήσουμε τη μεταβολή της περιοχής εξάπλωσης της στο μέλλον:

##==================================
## Load the the data
##==================================

## Current
Petromarula_pinnata_bin_proj_current <- stack('Petromarula.pinnata/proj_current/proj_current_Petromarula.pinnata_ensemble_TSSbin.img')
names(Petromarula_pinnata_bin_proj_current) <- c('ca', 'm', 'wm')

## 2070
Petromarula_pinnata_bin_proj_2070_CC85 <- stack('Petromarula.pinnata/proj_2070_CC85/proj_2070_CC85_Petromarula.pinnata_ensemble_TSSbin.img')
names(Petromarula_pinnata_bin_proj_2070_CC85) <- c('ca', 'm', 'wm')

##=================================
## Species Range Change
##=================================

SRC_current_2070_CC85 <- BIOMOD_RangeSize( Petromarula_pinnata_bin_proj_current,
                                           Petromarula_pinnata_bin_proj_2070_CC85 )
SRC_current_2070_CC85$Compt.By.Models

##    Loss Stable0 Stable1 Gain PercLoss PercGain SpeciesRangeChange
## ca  103     303      27   33    79.23    25.39             -53.85
## m   106     330      12   18    89.83    15.25             -74.58
## wm  100     297      31   38    76.34    29.01             -47.33
##    CurrentRangeSize FutureRangeSize.NoDisp FutureRangeSize.FullDisp
## ca              130                     27                       60
## m               118                     12                       30
## wm              131                     31                       69

Petromarula_pinnata_src_map <- SRC_current_2070_CC85$Diff.By.Pixel
names(Petromarula_pinnata_src_map) <- c("ca cur-2070", "m cur-2070", "wm cur-2070")

##=================================
## Vizualisation
##=================================

my.at <- seq(-2.5,1.5,1)

myColorkey <- list(
                   ## color change
                   at = my.at, 
                   labels = list(
                     
                     ## labels
                     labels = c("lost", "pres", "abs","gain"),
                     
                     ## legend
                     at = my.at[-1]-0.5 
                   ))

rasterVis::levelplot( SRC_current_2070_CC85$Diff.By.Pixel,
                      main = "Petromarula pinnata range change",
                      colorkey = myColorkey,
                      layout = c(1,3),
                      )

6 Εργασία για το σπίτι

Έχετε διορία δεκατεσσάρων (14) ημερών να στείλετε την εργασία σας σε μορφή PDF²³ σε αυτό το e-mail.
1. Κατεβάστε γεωγραφικά δεδομένα για την Ελλάδα.
2. Δημιουργήστε ξεχωριστά αρχεία για την Κρήτη, την Αττική, τον Άθω, την Μακεδονία και τη Θεσσαλία.
3. Κατεβάστε τα απαραίτητα κλιματικά δεδομένα.
4. Κατεβάστε υψομετρικά δεδομένα για την Ελλάδα.
5. Υπολογίστε τις βασικές στατιστικές παραμέτρους για όλα τα αρχεία τα οποία δημιουργήσατε.
6. Υπολογίστε τον συντελεστή συσχέτισης για όλα τα περιβαλλοντικά δεδομένα για όλα τα αρχεία που έχετε δημιουργήσει. Εάν υπάρχουν κάποια ζεύγη ανεξάρτητων μεταβλητών με σ.σ. \(>0.7\), θα τις χρησιμοποιήσετε όλες;
7. Κατεβάστε δεδομένα θέσης για τα είδη Asperula pubescens, Centaurea idaea και Scutellaria sieberi.
8. Παραθέστε πληροφορίες σχετικά με τα είδη αυτά (Καθεστώς κινδύνου/προστασίας, καθεστώς ενδημισμού, ταξινομικές πληροφορίες). Προχωρήστε στις απαραίτητες ενέργειες, έτσι ώστε να διαπιστώσετε εάν τα δεδομένα αυτά εντοπίζονται στη σωστή γεωγραφικά περιοχή.
9. Πόσες θέσεις εμφάνισης αριθμούσαν τα είδη αυτά πριν και μετά την διαδικασία αραίωσης;
10. Φορτώστε τα τελικά δεδομένα θέσης για κάθε είδος για το οποίο έχετε κατεβάσει δεδομένα θέσης.
11. Αλλάξτε τα rasterStack αρχεία με αυτά για την Κρήτη σε όλα τα απαραίτητα σημεία του κώδικα. Αλλάξτε επίσης το όνομα του είδους σας. Αλλάξτε τον αριθμό των σετ σε 10 και τον αριθμό των ψευδοαπουσιών²⁴. Αλλάξτε τον αριθμο των φορών που θα τρέξει η ανάλυση σε 10.
12. Ποιά είναι η σημαντικότερη μεταβλητή σε κάθε αλγόριθμο; Συμπίπτουν μεταξύ των αλγορίθμων;
13. Ποιό είναι το καλύτερο μοντέλο σύμφωνα με τους δείκτες αξιολόγησης TSS και KAPPA;.
14. Προβάλετε τα μοντέλα σας στο μέλλον. Παραθέστε τους σχετικούς χάρτες.
15. Πώς μεταβάλλεται η περιοχή εξάπλωσης των ειδών αυτών στο μέλλον;²⁵ Τι συμβαίνει στην περίπτωση της μηδενικής και της πλήρους διασποράς, αντίστοιχα;
16. Παραθέστε το σχετικό γράφημα.

7 Appendix: R code

##===========================================================================##
## Do NOT run this code
##===========================================================================##
library(knitr)
library(rmdformats)


## Global options
options(max.print="75")
opts_chunk$set(echo=TRUE,
                 cache=TRUE,
               prompt=FALSE,
               tidy=TRUE,
               # comment=NA,
               message=FALSE,
               warning=FALSE,
               eval = TRUE)
opts_knit$set(width=75)
##===========================================================================##
pre {
  max-height: 400px;
  overflow-y: auto;
}

pre[class] {
  max-height: 400px;
}
##===========================================================================##
## Install the main packages
##===========================================================================##
install.packages(c('tidyverse', 'rgbif', 'raster', 'data.table',  
                   'rgeos', 'rgdal', 'gridExtra', 
                   'dismo', 'biomod2', 'sf', 'usdm',
                   'speciesgeocodeR', 'pacman', 'spThin',
                   'CoordinateCleaner', 'ggspatial', 'scales',
                   'scrubr', 'mapr', 'rasterVis', 'ade4',
                   'viridis', 'biogeo'), dependencies = T)
##===========================================================================##
##===========================================================================##
## Load them
##===========================================================================##
pacman::p_load(tidyverse, rgbif, raster, usdm,  
                   rgeos, rgdal, usdm, biomod2, gridExtra, 
                   countrycode, devtools, sf, spThin, biogeo,
                   speciesgeocodeR, pacman, dismo,
                   CoordinateCleaner, ggspatial, scales,
                   scrubr, mapr, rasterVis, data.table, ade4, viridis)
##===========================================================================##
##===========================================================================##
## Download the data for Greece
##===========================================================================##
## There are three levels available, you are free to explore them
Greece <- getData('GADM', country = 'GRC', level = 3)
##===========================================================================##
##===========================================================================##
## Plot Greece
##===========================================================================##
plot(Greece)
##===========================================================================##
## First we need to see what is inside the object we created
##===========================================================================##
Greece
names(Greece)
Crete <- subset(Greece, NAME_1 == "Crete")
##===========================================================================##
## Plot Crete
##===========================================================================##
plot(Crete)

spp_petromarula <- name_suggest(q = 'Petromarula pinnata', rank='species',limit = 10000)
spp_petromarula
##===========================================================================##
## Download data
##===========================================================================##
sp_occ <- occ_search(taxonKey = spp_petromarula$data$key[1], 
                     country ='GR',
                     hasCoordinate = T,
                     limit = 1000)
##===========================================================================##


##===========================================================================##
## To prevent any problem with the pathway it is a good practice to remove blank
## spaces from species names
##===========================================================================##
sp_occ$data$name <- sub(" ", "_", sp_occ$data$name)
##===========================================================================##


##===========================================================================##
## Total number of occurrences
##===========================================================================##
sort(table(sp_occ$data$name), decreasing = T)
##===========================================================================##
## Clean the coordinates
##===========================================================================##
flags_gbif <- clean_coordinates(x = sp_occ$data, 
                           lon = "decimalLongitude",
                           lat = "decimalLatitude", 
                           countries = "countryCode", 
                           species = "species",
                           tests = c("capitals", 
                                     "centroids",
                                     "equal", 
                                     "gbif", 
                                     "zeros",
                                     "seas"), 
                           seas_ref = Crete)
##===========================================================================##


##===========================================================================##
## Keep only those that are OK
##===========================================================================##
flags_gbif <- flags_gbif %>% 
  mutate(Data = ifelse(.summary == 'TRUE', 
                       'Clean',
                       'Flagged'),
         NAME_1 = 'Crete',
         Dataset = 'GBIF',
         longitude = decimalLongitude,
         latitude = decimalLatitude,
         scientific_name  = 'Petromarula pinnata') %>% 
  dplyr::filter(str_detect(Data, 'Cle')) %>% 
  dplyr::select(scientific_name,
                longitude,
                latitude,
                Dataset) 
##===========================================================================##


##===========================================================================##
## Create a duplicate for safekeeping
##===========================================================================##
flags_gbif_dpl <- flags_gbif
##===========================================================================##
dir.create("./Species", recursive = TRUE, showWarnings = FALSE)
dir.create("./Species/Petromarula pinnata", recursive = TRUE, showWarnings = FALSE)

saveRDS(flags_gbif, "./Species/ Petromarula pinnata.rds")
##============================================================================##
## First, let's download current climate data for the whole planet
##============================================================================##
current_climate <- getData('worldclim', var = 'bio', res = 2.5)
##============================================================================##


##============================================================================##
## Then, let's download future climate data
##============================================================================##
future_climate <- getData('CMIP5', var = 'bio', res = 2.5, 
                          rcp = 85, model = 'CC', year = 70)

##============================================================================##
names(current_climate) <- c('bio_1', 'bio_2', 'bio_3', 'bio_4', 'bio_5', 
                               'bio_6', 'bio_7', 'bio_8', 'bio_9', 'bio_10', 
                               'bio_11', 'bio_12', 'bio_13', 'bio_14', 'bio_15', 
                               'bio_16', 'bio_17', 'bio_18', 'bio_19')

names(future_climate) <- c('bio_1', 'bio_2', 'bio_3', 'bio_4', 'bio_5', 
                               'bio_6', 'bio_7', 'bio_8', 'bio_9', 'bio_10', 
                               'bio_11', 'bio_12', 'bio_13', 'bio_14', 'bio_15', 
                               'bio_16', 'bio_17', 'bio_18', 'bio_19')
gain(current_climate$bio_1) = 0.1

gain(future_climate$bio_1) = 0.1
altitude <- getData('alt', country = 'GRC', mask = TRUE) 
##============================================================================##
## First for Greece
##============================================================================##
current_climate_crop <- crop(current_climate, 
                            Greece, 
                            snap = 'in') %>% 
  mask(., Greece)

altitude <- crop(altitude, 
                 Greece, 
                 snap = 'in') %>% 
  mask(., Greece) %>% 
  
  ## We are forcing the extent of altitude to be exactly the same as that of the
  ## climate data
  resample(., current_climate_crop, 'bilinear')  
  
future_climate <- crop(future_climate, 
                            Greece, 
                            snap = 'in') %>% 
  mask(., Greece)  
  
##============================================================================##


##============================================================================##
## Merge (stack) climate and altitude
##============================================================================##
current_climate_crop <- stack(current_climate_crop, altitude)
future_climate_crop <- stack(future_climate, altitude)
##============================================================================##


##============================================================================##
## Then for Crete
##============================================================================##
current_climate_Crete <- crop(current_climate_crop,
                              Crete,
                              snap = 'in') %>% 
  mask(., Crete)

future_climate_Crete <- crop(future_climate_crop,
                              Crete,
                              snap = 'in') %>% 
  mask(., Crete)

##============================================================================##
gplot(current_climate_Crete[[1]]) +
  geom_raster(aes(fill = value)) + 
  scale_fill_gradientn(colours = c("brown",
                                   "red", 
                                   "yellow", 
                                   "darkgreen", 
                                   "green")) + 
  coord_equal() + 
  xlab("Longitude") +
  ylab("Latitude") + 
  ggtitle("Temperature")

histogram(current_climate_Crete[[1]])

levelplot(current_climate_Crete[[1]])
cellStats(current_climate_crop, "quantile")
cellStats(current_climate_crop, "mean")
cellStats(current_climate_crop, "sd")
subset_GR <- current_climate_crop[[1]] < 15 & current_climate_crop[[1]] > 10.1

plot(subset_GR)
altitude
current_climate_crop
agg_altitude <- aggregate(altitude, 
                 5, ## Why choose 5 and not any other number?
                 fun = mean) 
altitude_rsmp <- resample(agg_altitude, current_climate_crop, method = "bilinear")
## ============================================================================##
## Merge (stack) climate and altitude
## ============================================================================##
current_climate_crop <- stack(current_climate_crop, altitude_rsmp)
future_climate_crop <- stack(future_climate_crop, altitude_rsmp)
## ============================================================================##
## Creating a folder to save our data
dir.create("./Data/", recursive = TRUE, showWarnings = FALSE)

## Save our data in rds format
saveRDS(current_climate_crop, "./Data/Environmental variables and Elevation.rds")
saveRDS(future_climate_crop, "./Data/Environmental variables and Elevation 2070.rds")

## Load our data from rds format
var_GR <- readRDS("./Data/Environmental variables and Elevation.rds")
var_GR_70 <- readRDS("./Data/Environmental variables and Elevation 2070.rds")
uncorrelated_vars <- vifcor(current_climate_Crete %>% as.data.frame(), th = 0.7)  
uncorrelated_vars <- uncorrelated_vars @results$Variables[uncorrelated_vars@results$VIF <= 5]

current_climate_Crete <- stack(subset(current_climate_Crete, 
                                      uncorrelated_vars))

future_climate_Crete <- stack(subset(future_climate_Crete, 
                                     uncorrelated_vars))

flags_gbif <- flags_gbif %>% as.data.frame()

coordinates(flags_gbif) <- c("longitude", "latitude")

pts.clim <- raster::extract(current_climate_Crete, flags_gbif, method = "bilinear")

sp_occ_clim <- data.frame(cbind(coordinates(flags_gbif), pts.clim, flags_gbif @data))
map.ext <- extent(Crete)

plot(current_climate_crop[[20]], col = terrain.colors(100), ext = map.ext)

plot(flags_gbif, pch = 16, cex = 0.5, add = T)

plot(Crete, add = T)
## Min. distance: 1000 m
## Number of reps: 1000
spoccs <- flags_gbif_dpl %>% 
  as.data.frame() %>%
  dplyr::rename(x = longitude,
         y = latitude,
         Species = scientific_name) %>% 
  addmainfields(species = "Species")
                                      
thinned <- thin(loc.data = spoccs, 
                lat.col = "y", 
                long.col = "x", 
                spec.col = "Species",
                thin.par = 5, 
                reps = 1000, 
                verbose = T, 
                out.dir = "./Thinned", 
                locs.thinned.list.return = TRUE,
                write.files = TRUE, out.base = "Petromarula pinnata thinned")
thinned_1 <- readxl::read_excel("G:/Academic/Lessons/EMB/Labs/Labs/Fourth_Lab/Fourth_Lab/Thinned/Petromarula pinnata thinned_thin1.xlsx")
sp_formated_data <- BIOMOD_FormatingData(resp.var = rep(1, nrow( thinned_1 ) ),
                                         
                                    ## Our rasterStack with the uncorrelated
                                    ## variables
                                    expl.var = current_climate_Crete,
                                    
                                    ## The names of our coords
                                    resp.xy = thinned_1[,c('x', 'y')],
                                    
                                    ## Species' name
                                    resp.name = "Petromarula pinnata",
                                    
                                    ## Number of pseudoabsences' sets
                                    PA.nb.rep = 2,
                                    
                                    ## Number equaling that of the presences
                                    PA.nb.absences = nrow(thinned_1), 
                                    
                                    ## You can change this to 'disk' if you'd
                                    ## like
                                    PA.strategy = 'random',
                                    
                                    ## Min distance from occurrences
                                    PA.dist.min  =  NULL,
                                    
                                    ## Max distance from occurrences
                                    PA.dist.max  =  NULL) 

sp_formated_data
plot(sp_formated_data)
modelling_options <- BIOMOD_ModelingOptions()
## Always change the name of this object, so as to match the name of your species
Petromarula_pinnata_models <- BIOMOD_Modeling( data = sp_formated_data,
                                       
                                       
                                       ## Here, we indicate which models we are
                                       ## going to run (remember that they need
                                       ## to be from the same group P/A ή P/O).
                                       ## Here, we are going to use some models
                                       ## belonging to the P/A group
                                       ## ('GLM','GBM','GAM','CTA','ANN','SRE',
                                       ## 'FDA','MARS','RF')
                                       models  =  c('RF', 'CTA'),
                                       
                                       models.options = modelling_options,
                                       
                                       ## How many times we will run the
                                       ## analysis in order to evaluate its'
                                       ## results [min: 10 - usually: 100]
                                       NbRunEval = 1,
                                       
                                       ## Training-Testing data split
                                       DataSplit = 80, 
                                       
                                       ## How many times we will run the
                                       ## analysis in order to evaluate the
                                       ## importance of each IV [min: 10]
                                       VarImport = 1,
                                       
                                       ## Evaluation criteria
                                       models.eval.meth  =  c('TSS','ROC','KAPPA'), 
                                       
                                       ## Whether or not we will save the results
                                       SaveObj  =  TRUE,
                                       
                                       ## Should the models be re-scaled in 0-1000 range?
                                       rescal.all.models  =  TRUE,
                                       
                                       ## Should our models be weighted based on all the data?
                                       do.full.models  =  TRUE, 
                                       
                                       ## The name of our analysis
                                       modeling.id  =  "allmodels") 

## Create some folders to store our results
dir.create('./BIOMOD2 Modelling Outpout/', recursive = TRUE, showWarnings = FALSE)

dir.create('./BIOMOD2 Modelling Outpout/Species/', recursive = TRUE, showWarnings = FALSE)

dir.create('./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata', recursive = TRUE, showWarnings = FALSE)

dir.create('./BIOMOD2 Modelling Outpout/Species/ Petromarula_pinnata/Current', recursive = TRUE, showWarnings = FALSE)


## Save the results in rds format
saveRDS(Petromarula_pinnata_models, './BIOMOD2 Modelling Outpout/Species/ Petromarula_pinnata/Current/ Petromarula_pinnata models.rds')
Petromarula_pinnata_scores <- get_evaluations(Petromarula_pinnata_models)


(Petromarula_pinnata_var_import <- get_variables_importance(Petromarula_pinnata_models))


options(digits = 2)
apply(Petromarula_pinnata_var_import, c(1,2), mean)

Petromarula_pinnata_ensemble_models <- BIOMOD_EnsembleModeling( modeling.output = Petromarula_pinnata_models,
                                                   em.by = 'all',
                                                   chosen.models  =  'all',
                                                   
                                                   # You can change that
                                                   eval.metric = 'TSS', 
                                                   
                                                   # Keep only the models with
                                                   # TSS >= 0.8
                                                   eval.metric.quality.threshold = 0.8, 
                                                   models.eval.meth = c('KAPPA','TSS','ROC'),
                                                   prob.mean  =  T, 
                                                   committee.averaging = TRUE,
                                                   prob.mean.weight  =  T, 
                                                   prob.mean.weight.decay  =  'proportional',
                                                   VarImport = 0 )




dir.create('./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Ensemble',
           recursive=TRUE, showWarnings=FALSE)

saveRDS(Petromarula_pinnata_ensemble_models,
        './BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Ensemble/Petromarula_pinnata ensemble.rds')
options(digits = 4)
(Petromarula_pinnata_ensemble_models_scores <- get_evaluations(Petromarula_pinnata_ensemble_models))

##===========================================================##
## Vizualisation 
##===========================================================##

g2 <- models_scores_graph( Petromarula_pinnata_ensemble_models, 
                            by = 'models',
                            metrics = c('TSS','KAPPA'))


Petromarula_pinnata_models_proj_current <- BIOMOD_Projection( modeling.output = Petromarula_pinnata_models,
                                                 new.env = current_climate_Crete,
                                                 proj.name = "current",
                                                 binary.meth = "TSS",
                                                 selected.models = 'all',
                                                 build.clamping.mask = TRUE,
                                                 output.format = ".img", 
                                                 do.stack = T )


Petromarula_pinnata_ensemble_models_proj_current <- 
  BIOMOD_EnsembleForecasting( EM.output = Petromarula_pinnata_ensemble_models,
                              projection.output = Petromarula_pinnata_models_proj_current,
                              binary.meth = "TSS",
                              output.format = ".img",
                              do.stack = T)


dir.create('./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/', 
           recursive=TRUE, showWarnings=FALSE)

dir.create('./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Individual', 
           recursive=TRUE, showWarnings=FALSE)

dir.create('./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Ensemble', 
           recursive=TRUE, showWarnings=FALSE)

saveRDS(Petromarula_pinnata_models_proj_current, 
        './BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Individual/Petromarula_pinnata Current Individual Projections.rds')

saveRDS(Petromarula_pinnata_ensemble_models_proj_current, 
        './BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Ensemble/Petromarula_pinnata Current Ensemble Projections.rds')


Petromarula_pinnata_proj_2070_cc85 <- BIOMOD_Projection( modeling.output = Petromarula_pinnata_models,
                                                   new.env = future_climate_Crete,
                                                   proj.name = "2070_CC85",
                                                   binary.meth = "TSS",
                                                   output.format = ".img",
                                                   do.stack = T )


dir.create('./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070', 
           recursive=TRUE, showWarnings=FALSE)

dir.create('./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/', 
           recursive=TRUE, showWarnings=FALSE)

dir.create('./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Individual', 
           recursive=TRUE, showWarnings=FALSE)

dir.create('./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Ensemble', 
           recursive=TRUE, showWarnings=FALSE)

saveRDS(Petromarula_pinnata_proj_2070_cc85, 
        './BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Individual/Petromarula_pinnata 2070 Individual Projections.rds')

Petromarula_pinnata_ensemble_models_proj_2070_cc85 <- 
  BIOMOD_EnsembleForecasting( EM.output = Petromarula_pinnata_ensemble_models,
                              projection.output = Petromarula_pinnata_proj_2070_cc85,
                              binary.meth = "TSS",
                              output.format = ".img",
                              do.stack = T )
saveRDS(Petromarula_pinnata_ensemble_models_proj_2070_cc85,
      './BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Ensemble/Petromarula_pinnata 2070 Ensemble Projections.rds')
Petromarula_pinnata_ensemble_models_proj_2070_cc85 <- readRDS('./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Ensemble/Petromarula_pinnata 2070 Ensemble Projections.rds')
stk_Petromarula_pinnata_ef_2070_cc85 <- get_predictions(Petromarula_pinnata_ensemble_models_proj_2070_cc85)

stk_Petromarula_pinnata_ef_2070_cc85 <- subset(stk_Petromarula_pinnata_ef_2070_cc85,
                                  grep('EMmean|EMwmean',
                                       names(stk_Petromarula_pinnata_ef_2070_cc85)))

names(stk_Petromarula_pinnata_ef_2070_cc85) <- sapply(strsplit(names(stk_Petromarula_pinnata_ef_2070_cc85),
                                                  '_'),
                                         getElement, 2)

Petromarula_pinnata_coords <- thinned_1 %>% dplyr::select(x, y)
library(viridis)
levelplot(stk_Petromarula_pinnata_ef_2070_cc85,
          main = 'Petromarula pinnata ensemble projections\nin 2070 with cc85',
          col.regions = viridis(256)) +
  layer(sp.points(SpatialPoints(Petromarula_pinnata_coords), pch=20, col=1))
##==================================
## Load the the data
##==================================

## Current
Petromarula_pinnata_bin_proj_current <- stack('Petromarula.pinnata/proj_current/proj_current_Petromarula.pinnata_ensemble_TSSbin.img')
names(Petromarula_pinnata_bin_proj_current) <- c('ca', 'm', 'wm')

## 2070
Petromarula_pinnata_bin_proj_2070_CC85 <- stack('Petromarula.pinnata/proj_2070_CC85/proj_2070_CC85_Petromarula.pinnata_ensemble_TSSbin.img')
names(Petromarula_pinnata_bin_proj_2070_CC85) <- c('ca', 'm', 'wm')

##=================================
## Species Range Change
##=================================

SRC_current_2070_CC85 <- BIOMOD_RangeSize( Petromarula_pinnata_bin_proj_current,
                                           Petromarula_pinnata_bin_proj_2070_CC85 )
SRC_current_2070_CC85$Compt.By.Models

Petromarula_pinnata_src_map <- SRC_current_2070_CC85$Diff.By.Pixel
names(Petromarula_pinnata_src_map) <- c("ca cur-2070", "m cur-2070", "wm cur-2070")

##=================================
## Vizualisation
##=================================

my.at <- seq(-2.5,1.5,1)

myColorkey <- list(
                   ## color change
                   at = my.at, 
                   labels = list(
                     
                     ## labels
                     labels = c("lost", "pres", "abs","gain"),
                     
                     ## legend
                     at = my.at[-1]-0.5 
                   ))

rasterVis::levelplot( SRC_current_2070_CC85$Diff.By.Pixel,
                      main = "Petromarula pinnata range change",
                      colorkey = myColorkey,
                      layout = c(1,3),
                      )
##===========================================================================##
## Do NOT run this code
##===========================================================================##
labs = knitr::all_labels(echo = TRUE)
##===========================================================================##
klippy::klippy()

8 Βιβλιογραφία

Araújo, M.B. & New, M. (2007) Ensemble forecasting of species distributions. Trends Ecol. Evol., 22, 42–47.

Guisan, A., Thuiller, W., & Zimmermann, N.E. (2017) Habitat Suitability and Distribution Models with Applications in R. Cambridge University Press,

McSweeney, C.F., Jones, R.G., Lee, R.W., & Rowell, D.P. (2015) Selecting CMIP5 GCMs for downscaling over multiple regions. Clim. Dyn., 44, 3237–3260.

Υπάρχουν δύο τρόποι για να το κάνετε αυτό γρήγορα:
1. library(easypackages) packages('package_name_1', 'package_name_2')
2. install.packages(c('package_name_1', 'package_name_2'), dependecies = TRUE)↩︎
Ή με αυτήν την εντολή εάν θέλουμε μόνο την Ελλάδα: ISO_countries <- getData("ISO3") %>% as.data.frame %>% filter(NAME=="Greece")↩︎
Ένα τυπικό παράδειγμα είναι ο Άτλαντας της Χλωρίδας του Αιγαίου ή ο Άτλαντας της Χλωρίδας της Ευρώπης, όμως στην περίπτωση αυτή υπεισέρχονται άλλα προβλήματα τα οποία σχετίζονται με την ανάλυση των εν λόγω χαρτών (π.χ., ο Άτλαντας της Χλωρίδας της Ευρώπης είναι σε ανάλυση \(50 \times 50 km^2\))↩︎
Δεν είναι τα μοναδικά όμως: υπάρχουν πολλά άλλα, όπως το dismo και το spocc ↩︎
Κάθε ένα από αυτά τα πακέτα, έχει τουλάχιστον μια σύντομη περιγραφή (αγγλιστί: vignette), την οποία είναι ΑΠΑΡΑΙΤΗΤΟ να διαβάσετε - ειδικά αυτές του rgbif↩︎
Με την εντολή ?rgbif::name_suggest() μπορείτε να δείτε τα arguments που δέχεται η εντολή αυτή και τι μπορείτε να αλλάξετε - σας συνιστώ εντόνως να το κάνετε, αλλά και για κάθε εντολή που χρησιμοποιείτε↩︎
Καλό είναι να τις διαβάσετε, γιατί θα σας βοηθήσει και στην εργασία σας↩︎
Και εγώ στη θέση σας, το ίδιο θα έκανα↩︎
Γιατί;↩︎
Γιατί όμως να μην αλλάξουμε την ανάλυση των περιβαλλοντικών δεδομένων;↩︎
Αυτό είναι το λεγόμενο hard boundary - άλλοι ερευνητές εφαρμόζουν το όριο \(> 0.8\) ή ακόμα και \(> 0.9\), αλλά είναι επιστημονικά ασφαλέστερο να εφαρμόζει κανείς το αυστηρότερο όριο↩︎
Εδώ εφαρμόζουμε το αυστηρότερο όριο - μπορείτε να το αλλάξετε, εάν επιθυμείτε↩︎
Σε απόσταση \(<1000m\) όταν εργαζόμαστε με χάρτες υποβάθρου κλίμακας 1 x 1 km² (ή ανάλυσης 30 arc-seconds)↩︎
Θα σας βοηθήσει το vignette του εν λόγου πακέτου↩︎
Προτρέξτε στους Guisan et al. (2017) για περισσότερες λεπτομέρειες - εάν χρειάζεστε περισσότερη βιβλιογραφία, απευθυνθείτε σε εμένα (εάν φυσικά το επιθυμείτε)↩︎
Όπως θα διαπιστώσετε παρακάτω, θα εφαρμόσουμε μια τεχνική η οποία ονομάζεται ensemble modelling. Σύμφωνα με αυτή, είναι προτιμότερο και επιστημονικά ορθότερο να μην χρησιμοποιεί κανείς ένα μόνο μοντέλο (π.χ., GLM ή τo MaxEnt)↩︎
Προσέξτε ότι πρόκειται για αρχείο xlsx και όχι csv - αυτό σημαίνει ότι άνοιξα στο excel το αρχείο csv, χρησιμοποίησα την εντολή data to columns και το έσωσα ως αρχείο xlsx↩︎
Δεν έχει σημασία ποιό θα διαλέξουμε - έχουν όλα την ίδια πιθανότητα↩︎
Θα μπορούσε κάποιος να αντιτείνει ότι μπορούμε να χρησιμοποιήσουμε το πακέτο ecospat όταν \[\frac{\sum{Ανεξάρτητες,\, μη\, συγγραμικές\, μεταβλητές}}{\sum{Θέσεις\, εμφάνισης}} \geq 0.1\]↩︎
Στο βιβλίο των Guisan et al. (2017) αναφέρεται λανθασμένα ότι ο αριθμός των pseudoabsences πρέπει να είναι τουλάχιστον 10000 - ο αριθμός αυτός έχει να κάνει με το MaxEnt, το οποίο ανήκει στην ομάδα των P/O μοντέλων↩︎
Καλό είναι να είναι τουλάχιστον 10 τα set αυτά↩︎
Θυμηθείτε να με ρωτήσετε τι σημαίνουν οι υπόλοιπες μεταβλητές↩︎
Διαφορετικά δεν θα γίνει δεκτή και δεν θα βαθμολογηθείτε↩︎
Πάντα σύμφωνα με το είδος σας↩︎
Παραθέστε ΟΛΑ τα στοιχεία από τον σχετικό πίνακα↩︎

Delivered to you by Dr. Kostas Kougioumoutzis

kkougiou@aua.gr